/***************************************************************************
*   Copyright (C) 2009 by Steven Waslander   *
*   stevenw@uwaterloo.ca   *
*                                                                         *
***************************************************************************/


#include <iostream>
#include "sampleCode/Mathutils.h"

using namespace std;

void DebugTest1(void);
void DebugTest2(void);
void DebugTest3(void);

int main()
{

  DebugTest1();
  DebugTest2();
  DebugTest3();
  return 0; 
}


void DebugTest3(void)
{
  Vec<double> x(4);
  x(0)=1;
  x(1)=2;
  x(2)=3;
  x(3)=4;
  Vec<double> y(x);
  cout << "Testing Vector library" << endl;
  cout << "Define vector x: " << endl << x << endl;
  cout << "Copy x to y: " << endl << y << endl;
  cout << "Sum x+y: " << endl << x+y << endl;
  cout << "Difference x-y: " << endl << x-y << endl;
  cout << "Product x'*y: " << endl << x.Transpose()*y << endl;
}

void DebugTest2(void)
{
  int N = 4;
  int M = 4; 
  int d;
  Matrix<int> indx(4,1);
  Matrix<double> b(4,1);
  b(0,0) = 0.0;
  b(1,0) = 0.0;
  b(2,0) = 1.0;
  b.Randn(0.1,3.0);
  cout << "Testing Matrix library" << endl;

  Matrix<double> A(N,M);
  for (int i=0;i<N;++i) {
    for (int j=0;j<M;++j) {
      A(i,j) = (i*i-j*j*j)/(i+1)+3;
    }
  }

  cout << "Define a matrix, A:" << endl << A << endl;
  cout << "Transpose of A, A':" << endl << A.Transpose() << endl;
  cout << "Trace(A): " << A.Trace() << endl;
  seedRand();
  A.LuDecomp(indx,&d);
  cout << "LU Decomposition of current A:" << endl << A << endl;
  cout << "LU Decomposition Index:" << endl << indx << endl;
  cout << "Define a vector b:" << endl << b << endl;
  A.LuSolve(indx, b);
  cout << "Solve Ax=b with LU Decomposition:" << endl << b;
  cout << "Random A, Uniform:" << endl <<   A.Rand(1.0,2.0) << endl;
  cout << "Random A, Gaussian:" << endl << A.Randn(0.1,2.0) << endl;
  A.Inverse();
  cout << "A inverse:" << A;  
  cout << "Det(A):" << endl << A.Det() << endl;
}

void DebugTest1(void)
{
  int N = 3;
  int M = 3;
  cout << "Testing Matrix library" << endl;

  Matrix<double> A(N,M);
  for (int i=0;i<N;++i) {
    for (int j=0;j<M;++j) {
      A(i,j) = (i+j)/2.0;
    }
  }
  cout << "A Matrix" << endl;
  cout << A;

  Matrix<double> B(M,1);
  for (int i=0;i<M;++i) {
    B(i,0) = 0.3 - i*1.1;
  }
  cout << "B Matrix" << endl;
  cout << B;

  Matrix<double> C(A);
  cout << "C, Copy of A" << endl;
  cout << C;
  cout << "Matrix Addition: A+=C" << endl;
  A += C;
  cout << A;
  cout << "Matrix Subtraction C-=A" << endl;
  C -=A;
  cout << C;
  cout << "Matrix Multiplication A*B" << endl;
  cout << A*B;
  
  ofstream outmat;
  outmat.open("outmat.txt");
  if (!outmat)
    cout << "File open failed" << endl;

  outmat << A;
  outmat.close();

  ifstream inmat;
  inmat.open("outmat.txt");
  if (!inmat)
    cout << "File open failed" << endl;

  Matrix<double> D(N,M);
  inmat >> D;
  cout << "D: read from outmat.txt" << endl;
  cout << D;
   
  inmat.close();

}

